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It has been shown previously that a finite-length topological insulator nanowire, proximity-coupled 
to an ordinary bulk s-wave superconductor and subject to a longitudinal applied magnetic field, 
realizes a one-dimensional topological superconductor with an unpaired Majorana fermion (MF) 
localized at each end of the nanowire. Here, we study the stability of these MFs with respect to var- 
ious perturbations that are likely to occur in a physical realization of the proposed device. We show 
that the unpaired Majorana fermions persist in this system for any value of the chemical potential 
inside the bulk band gap of order 300 meV in Bi2Se3 by computing the Majorana number. From 
this calculation, we also show that the unpaired Majorana fermions persist when the magnetic flux 
through the nanowire cross-section deviates significantly from half fiux quantum. Lastly, we demon- 
strate that the unpaired Majorana fermions persist in strongly disordered wires with fluctuations in 
the on-site potential ranging in magnitude up to several times the size of the bulk band gap. These 
results suggest this solid-state system should exhibit unpaired Majorana fermions under accessible 
conditions likely important for experimental study or future applications. 



I. INTRODUCTION 

In 1937, Ettore Majorana first shovifed that the com- 
plex Dirac equation can be separated into a pair of 
real wave equations, each of which is satisfied by real 
fermionic fields-. Such a real fermionic field, denoted by 
^, satisfies the property that ^ = A particle created 
by this field, known as a Majorana fermion, is therefore 
distinguished by the fact that it is its own antiparticld^^'^. 
Having many properties which make them interesting 
from the standpoint of fundamental science, while also 
being a possible platform for fault-tolerant, scalable 
quantum computation^ Majorana fermions are of 
tremendous interest to the condensed matter commu- 
nity. After intense effort, some proposa ls to realize Ma- 
jorana zero-modes made in recent year^^^^^ seem to be 
bearing fruit, with signatures of Majorana fermions al- 
ready being reportecpSHiZl Qf the many devices proposed 
for harbouring Majorana fermion^^, however, virtually 
all face considerable experimental challenges in achiev- 
ing the conditions necessary for Majorana fermion emer- 
gence. Additional hurdles are associated with the control 
and manipulation of MFs which is necessary for harness- 
ing their potential for quantum computation. Thus, even 
if Majorana fermions have indeed been conclusively ob- 
served, there remains a need for more accessible plat- 
forms with which to realize MFs sufficiently robust for 
applications. 

The purpose of this paper is to present results on 
the stability of MFs in a solid-state device previously 
predictecpSl to host these quasiparticle excitations. The 
device, depicted schematically in Fig. [TJ consists of a 
nanowire fashioned out of a strong topological insulator 
(STI), such as Bi2Se3 or Bi2Te2Se, placed in contact with 
an ordinary s-wave superconductor (SC), subject to an 
applied magnetic field along the axis of the nanowire. 
We show that MFs are remarkably stable in this de- 
vice, making it unique amongst the many proposals for 
observing MFs in solid-state systems and a significant 
advancement towards study of MFs and development of 




FIG. 1: Schematic of the proposed device. Magnetic field B 
is applied along the axis of the wire taken to coincide with 
the 2-direction. 



MF-based technology. 

We note that Bi2Se3 nanowires and nanoribbons have 
been synthesized and can exhibit diverse morphologies 
controllable by growth conditionJi^l, Aharonov-Bohm 
(AB) oscillations in the longitudinal magneto-resistance 
of Bi2Se3 nanoribbons have also been observed, proving 
the existence of a coherent surface conducting channepSl. 
Studies of magneto-resistance of Bi2Se3 nanoribbons un- 
der a variety of magnetic field orientations also reveal a 
linear magneto-resistance that persists to room tempera- 
ture and is consistent with transport through topological 
surface state^^. Lastly, the superconducting proximity 
effect and possible evidence for Pearl vortices has been 
observed in Bi2Se3 nanoribbon^^H This experimental 
progress suggests our proposed device may be realized 
experimentally with relative ease. 

Stability of MFs in our proposed device is confirmed by 
showing that the degenerate quasiparticle ground state is 
separated from excited states by an energy gap close to 
the superconducting (SC) gap which can be as large as 
^ 10 meV, through study of a low-energy analytical the- 
ory and numerical study of a lattice-model Hamiltonian. 
We also compute the topological phase diagram for the 
system numerically to show that MFs exist in the system 
for any value of the chemical potential in the bulk band 
gap of the TI (for Bi2Se3, ^ 300 meV). Furthermore, we 
find that the topological phase corresponding to the pres- 
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ence of MFs persists even when the chemical potential is 
in the bulk conduction band, although, since we observe 
rapid collapse of the excitation gap in this regime, this 
result is of limited experimental relevance. 

These results also support additional explicit numeri- 
cal studies of the robustness of MFs against non-magnetic 
disorder also discussed, which show that the MFs persist 
in the presence of fluctuations in the on-site chemical 
potential in an explicit lattice model Hamiltonian on the 
order of the bulk band gap (300 meV). As such, previ- 
ous expectations that MFs would be robust against non- 
magnetic disorder according to Anderson's theorem, be- 
cause time reversal symmetry (TRS) holds in this device 
under operating conditional^], are here confirmed. 

In the most promising of other proposals, the MFs are 
protected by an SC gap of at most 1 meV, and the chemi- 
cal potential must also be tuned to lie within a window of 
the same siz^^^^. Although these requirements are pos- 
sible to achieve in experiments on individual wire^', such 
fine tuning will be difficult to replicate in more complex 
setups, i.e. those containing wire networks necessary for 
MF manipulatiorP31. Furthermore, other proposals are 
not predicted to possess MFs under TR-invariant con- 
ditions, meaning these devices are not expected to be 
robust against non-magnetic disorder. Therefore MFs 
constructed in these proposals might be too delicate to 
be useful in practical applications. Our results on the 
remarkable stability of MFs in the TI nanowire-based 
proposal therefore outline a practical route towards ap- 
plications based on the physics of Majorana fermions. 



II. MAJORANA FERMIONS 
A. Significance and relevance of Majorana fermions 

In the more than 70 years since Majorana's seminal 
papei^, it is only very recently that any experimental 
evidence of Majorana fermions has been obtained. Re- 
cent developments in topological states of matter have 
potentially already led to successful const ruction of Ma- 
jorana fermions in solid-state system^^, as quasiparti- 
cle excitationJ^. The race to experimentally study these 
quasiparticles has been especially heated because of the 
observation by Kitae\i2l! that Majorana fermions could 
be used as a platform for robust quantum computing. In 
this scheme, a quantum bit is stored in a Dirac fermion 
that has been teased apart into two Majorana fermions. 
If these two Majorana fermions are separated spatially 
from one another, then, whether this shared fermion is 
occupied or empty, it is distributed non-locally, and no lo- 
cal perturbation can measure this shared quantum biiP^. 

Furthermore, a system of N spatially-separated Majo- 
rana fermions is predicted to satisfy non-Abelian statis- 
tics, implying that such a system has an iV-quasiparticle 
ground state that is degenerate. This degeneracy allows 
adiabatic interchange of the quasiparticles, or braiding, 
to correspond to unitary operations on the ground state. 



For Majorana fermions, it has also been shown that the 
only way to perform unitary operations on the ground 
state - which could be used for computing - is by braiding, 
and these operations are dependent only on the topol- 
ogy of the braid. Since the system is in a topological 
phase when Majorana fermions are present, this degen- 
erate ground state is also separated from the rest of the 
spectrum by an energy gap known as the "minigap". If 
the temperature is much lower than the minigap, and 
the system is weakly perturbed using frequencies much 
smaller than the gap, the system evolves only within the 
ground state subspac^^H, 

All of these features combined mean that a system of 
spatially-separated Majorana fermions could be used as a 
quantum computer that is immune to the tremendous ob- 
stacle faced by most other proposed platforms for quan- 
tum computing known as decoherenc^^H Experimental 
confirmation of the existence of Majorana fermions is a 
crucial first step towards practical quantum computing, 
but it is imperative that platforms possessing robust Ma- 
jorana fermions under stable conditions be identified and 
developed. 



B. Other existing proposals for realizing Majorana 
fermions experimentally 

There is no shortage of proposals for realizing Ma- 
jorana fermions experimentally. Earlier suggestions for 
physical systems that support Majorana fermion states 
include fractional quantum Hall states at filling ly = 
and Helium-^PS. These ground-breaking propos- 
als are thought to be extremely challenging to realize 
experimentallj^, however. We will discuss the many 
other proposals and comment on the experimental chal- 
lenges they face below. 

2D topological insulators have long been proposed 
as platforms for realizing Majorana zero-modes, for in- 
stance, having the advantages of greatly facilitating 
Josephson-based Majorana detection, long considered to 
be smoking gun confirmation of the presence of Majorana 
zero- modes'*, as well as being unaffecte d by n on-magnetic 
disorder due to time-reversal invarianc j^^ l I and, in prin- 
ciple, possessing a la rge p airing gap exhibited by the 
parent superconductop2212U However, of many materials 
predicted to be 2D topological insulators^^, only one, 
HgTe, has been confirmed experimentally thus fai'^^ '^2', 
although there has also been some evidence recently that 
InAs/GaSb quant um w ells may also exhibit a topologi- 
cal insulator phasd^^^. 2D semiconductor heterostruc- 
tures have also shown promise as platforms for realizing 
Majorana zero-m odes, but face challenges due to small 
spin-orbit energie^^^Ell^ a need for difficult-to-engineer, 
high-quality interfaces, and limited tunabilit}^. 

An innovative proposal for realizing Majorana zero- 
modes in three dimensional topological insulators due to 
Fu and Kane exists^^j^ but this proposal, while ground- 
breaking, faces considerable challenges given that time- 
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reversal symmetry must be broken to achieve Majo- 
rana zero-modes, making the device vulnerable to non- 
magnetic disordeil^. There have also been many propos- 
als based on SU2RUO4, but even in the simplest of these 
proposals, the minigap protecting Majorana zero-modes 
from excited states is in the milliKelvin rang4^, and a 
beautiful proposal for realizing Kitaev's ID toy model 
along an ordinary |^ vortex line threading a layered spin- 
ful p + ip superconduct or lik ely to be SU2RUO4 currently 
faces the same problerrPEH. 

There is great interest in realizing Majorana zero- 
modes in one dimensional systems, because they have 
generally been predicted to remain separated from ex- 
cited states by a larger energy gap than in other 
proposal J^. Conventional ID wires with sizeable spin- 
orbit coupling, proximate to s-wave super conductors, 
and subject to modest mag netic fieldPE3 are seen as 
very promising platforms for first experimental realiza- 
tion of Majorana zero-mode^. These proposals must 
overcome numerous issues, however, such as position- 
ing of the chemical potential in a rather small inter- 
val of roughly 1 K over distances long compared to the 
wire's coherence lengtlP. This constraint could be re- 
laxed by applying larger magnetic fields, but this intro- 
duces other difhcultie^. Tuning of the chemical poten- 
tial could likely be even more difficult due to disorder- 
induced fluctuations in the chemical potential, since the 
topological phase corresponding to the presence of Ma- 
jorana zero-modes appears only at finite magnetic fields 
in these devices, so Anderson's theorem does not protect 
the gap against non-magnetic disorder, which i s always 
pair-breaking according to many previous studieP^'^SMH 
Further, since the ratio of Zeeman energy to spin-orbit 
energy is small for both wires made of InAs and InSlP^, 
disorder is likely to play a non-trivial rol^. Although 
there have been efforts to ameli orate this issue by elimi- 
nating an appli ed magnetic fielcP^^^ from the device or 
reducing itP^ES^ these approaches can also lead to com- 
plications that can poten tially cause the Majorana zero- 
modes to disappeapSSEZI, 

The above conventional ID wire proposals further face 
the challenge that multiple sub-bands are usually oc- 
cupied in these wires and gating into the lowest sub- 
band regime is potentially non-trivial, especially if these 
wires are in close proximity to a superconductor as 
proposecP. Multichannel wires have been shown to sup- 
port the ID topological superconductor state leading to 
Maj orana zer o-modes away from the lowest sub-band 
limiil^215 2 | 59|60| ^ but these systems still require some de- 
gree of gating, leading to proposals of increasing com- 
plexity involving regular arrays of superconducting is- 
lands in contact with the wir^^^^. Such work has even 
led to the ingenious proposal of a chain of quantum dots 
that would be bridged by superconducting islands®^, but 
reaching the regime where only a comparatively small 
number of quantum dots would be needed would re- 
quire very fine-tuning that would likely suffer from strong 
randomnesJ^. Carbon nanotubes, also suggested as hosts 



for Majorana zero-modes, face considerable challenges in 
reaching the spinless regime with proximity-induced pair- 
ing requirecP^t^, while proposals involving half-metallic 
ferromagnetic wire also face challenges, such as the need 
to couple to non-c entro symmetric superconductors with 
spin-orbit couplin^sZMll 

Despite these challenges, there have been promis- 
ing experimental results for Majorana fermions based 
on some of these proposals. Josephson effects at the 
surface of a variety of 3D topological insulators with 
superconducting electrodes have been observed^^ 69 74]^ 
While these experiments, and related Andreev conduc- 
tance measurementa^SHlIl show interesting and unusual 
features, these cannot be readily attributed to the sin- 
gle Majorana zero-mode (typically only one out of 10^ 
modes)^. 

The nanowire-based proposal of Lutchyn et aL^i^E^ and 
Oreg et al.^ has also led to convincing evidence for a 
Majorana zero-mode in an InSb nanowire as reported by 
Kouwenhoven and his group^^. Since then, theoretical 
worlP^ by Patrick Lee and collaborators has indicated 
that, under conditions for semiconducting wires with 
modest amounts of disorder relevant to Kouwcnhovcn's 
work, Majorana end-states are destroyed and do not give 
rise to quantized zero-bias peaks (ZBPs). At finite tem- 
peratures, furthermore, ZBPs of a non-topological origin 
are predicted to appear, leading to clusters of low-energy 
states localized near the wire end. These non-topological 
ZBPs are further anticipated to be typically stable with 
respect to variations of chemical potential and magnetic 
field, and appear and disappear under nearly identical 
conditions to those of true Majorana peaks. This work 
suggests caution is required in interpreting recent exper- 
iments to observe MZMs and that substantially longer 
and cleaner wires are required to conclusively observe 
MZMs. 

However, work by Tewari and Stanescu^ also indi- 
cates, for a smooth confinement potential at the ends 
of a semiconductor Majorana wire, emergence of zero 
bias conductance peaks corresponding to the topologi- 
cally trivial phase is necessarily accompanied by a sig- 
nature similar to closing of the bulk band gap. This 
gap closing signature in the end-of-wire local density of 
states was absent in the Kouwenh oven study, suggest- 
ing Kouwenhoven's group and otherd^^^ may have been 
successful in observing MZMs. If indeed Majorana zero- 
modes have finally be observed, however, there still re- 
mains a need for devices in which MFs can be realized un- 
der more accessible conditions, are robust, and can finally 
be manipulated for topologically-protected computation, 
motivating the results we present here on a proposal in 
which Majorana fermions occur under a wide-range of 
accessible conditions robustly. 
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III. TI NANOWIRE WITH MAGNETIC AND 
SUPERCONDUCTING ORDER 

A. Low-energy theory: normal state 

We begin by presenting the low-energy analytical the- 
ory of the devic^i^ in greater detail to facilitate later 
discussion of the novel results on stability, as this foun- 
dation is later used to understand the new results. 

First, we motivate the proposal with study of a cylin- 
drical TI nanowire proximity-coupled to a bulk s-wave 
SC as the greater symmetry of this system permits an- 
alytical study of the low-energy fermionic excitations on 
the surface of the nanowire. 

The low-energy fermionic excitations on the surface 
of the topological insulator are governed by the Dirac 
HamiltoniarpS 

ho ^ - [hW ■ h + n ■ {p X s) + {p X s) ■ h], (1) 

where n is a unit vector normal to the surface, p = — iV 
is the momentum operator and s is the vector of Pauli 
matrices in the spin space. We will also include the effect 
of a magnetic coating on the TI nanowire by adding an 
additional term, h^ = s • m, to the Hamiltonian. Later, 
we will show that this term is not necessary for Majorana 
zero-modes to emerge in the device, but its inclusion will 
be convenient in calculations. 

Let us now consider the specific case of a cylindrical 
topological insulator nanowire of radius R with magnetic 
field B applied along the z-axis as shown in Fig. [2] The 
unit vector n is then taken to be normal to the curved sur- 
face of the nanowire. To include a flux $ through the end 
of the wire (in the z direction) as proposed, we replace the 
momentum operator p in Eq. ([ij with tt = p — (e/c)A, 
where A — rj^oiz x r)/27rr^ is the vector potential and 
$0 is the flux quantum. Therefore, suppressing vh, we 




FIG. 2: Schematic of the device simplified for analytical study, 
in which a cylindrical TI nanowire is substituted for a more 
realistic TI nanowire with square cross-section. Magnetic field 
B is still applied along the axis of the wire taken to coincide 
with the «-direction. 



now have the Hamiltonian, 

= -'-I+(n X 7r)-s + s • m. (2) 
2r 

Taking m = mz, we can rewrite the Hamiltonian in cylin- 
drical coordinates as 

-^I+siksm{(j))-S2kcos{(l))-S3 (^^'^'A + +™'^3- 

(3) 

To diagonalize this Hamiltonian for an infinitely long 
wire, we exploit the translational and rotational sym- 
metries and write a solution tpki of the form 

V'fe(z,^) = e''^'e-''=^(^^,{^^^J (4) 
With this ansatz, our Hamiltonian is 

hki = S2k + S3[{l+^~Tj)/R + m]. (5) 
The spectrum Eki for m = 0, if vh is reinstated, is then 



Ek, = ±vh]^ fc2 + ^ . (6) 

Here k labels momentum eigenstates along the cylinder 
while I = 0,±1, ... is the angular momentum. We see 
that the spectrum has a gapless branch for r/ = n + |, 
where n is any integer {rj = ^/^o measures the magnetic 
flux through the wire cross section in the units of flux 
quantum $o — hc/e). The periodicity r] ^ rj + n with n 
integer in Eq. ^ reflects the expected <I>o-periodicity in 
the total flux. 

We now notice that, although the branches of Eki are 
doubly-degenerate for 77 = 0, the degeneracy is lifted for 
77 7^ 0. One can always find a value of the chemical po- 
tential /i that yields a single pair of non-degenerate Fermi 
points for 77 7^ 0, as illustrated in Fig. [3] or more gener- 
ally, an odd number of such pairs. According to the Ki- 
taev's criteriorP^ pairing induced by the proximity effect 
is then expected to drive the system into a topological 
phase. In the special case 77 = 1/2 the two lowest bands 
are non-degenerate, while all higher bands are doubly- 
degenerate, yielding an odd number of Fermi points at 
any /i in the bulk band gap. While the semiconductor 
wire proposa P'^ l ^^ l only possesses Majorana fermions for 
values of the chemical potential in a 1 meV interval, our 
proposal possesses Majorana fermions for any value of 
the chemical potential fi inside the 300 meV bulk band 
gap of Bi2Se3. Fine-tuning of the chemical potential is 
unnecessary in our device at 77 = 1/2 due to the specific 
pattern of degeneracies of the bands which is in turn pro- 
tected by the Kramers theorem. 

The surface Dirac Hamiltonian ([!]) is expected to be 
valid in the limit when the surface state penetration 
depth C is much smaller than the wire radius R. In the 
opposite limit of a thin wire, C > one could worry 
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FIG. 3: The normal state dispersion under conditions 
required for emergence of Majorana fermions for (a) 
t he Ra shba-coupled semiconductor quantum wire proposal 
^ ji3 | i4 |^ where the dispersion is shown without Zeeman cou- 
pling (dashed lines) and with Zeeman coupling, and (b) 
our topological insulator nanowire proposal, with doubly- 
degenerate bands shown as black and blue dashed lines. Green 
and pink horizontal lines represent the level of the chemical 
potential and a number to the right of a line indicates the 
number of Fermi points in the right half of the Brillouin zone 
at that value of the chemical potential. Vertical green lines 
indicate the interval inside which the chemical potential can 
be tuned to yield Majorana fermions in the corresponding SC 
state. 



that the wavefuncion overlap in the interior of the wire 
might lead to the formation of a gap, as happens e.g. in 
thin TI films. We study the limit of a thin wire in Ap- 
pendix A, based on a 3D effective model of a TI. The 
results of this study are interesting. We find that the 
gapless mode actually persists for an arbitrary radius R 
in the case when the magnetic flux $o/2 has the form of 
a (5-function centered at the axis of the cylinder. For a 
uniform magnetic flux a gap in the surface state opens up 
and its magnitude is proportional to {C,/ K)^ representing 
the amount of magnetic flux to which the wave function 
is exposed. Thus, in the case of a cylindrical wire, it is 
not the wavefunction overlap (which would lead to a gap 
g-R/C'j byl; ^hc amouut of T-breaking in the system 
that determines the gap. We note that our numerical 
simulations discussed in Sec. IV below indicate that for 
moderately thin wires, i.e. C being a significant fraction 
of R, the gapless state actually persists but now exists at 
somewhat higher magnetic flux. Also, the wires likely to 
be used in an experimentP^'^ are tens of nm thick and 
are thus in the thick- wire limit, ^ being typically just a 
few lattice spacings. 

Finally we note that the surface Dirac Hamiltonian 
^ represents the simplest possible model that neglects 
anisotropies present in real materials, such as Bi2Se3. 
Such anisotropies will necessarily lead to a spin texture 
that depends on the crystallographic direction of the sur- 
face which can be described by an effective 4x4 sur- 
face Hamiltonian.'^ It would be interesting to understand 
these effects for various wire geometries but we leave this 
problem for future study. We note that our lattice Hamil- 
tonian employed in Sec. IV includes the above mentioned 



anisotropies and the results based on it confirm all the 
essential features of our simple analytical model. 

B. Low-energy theory: Majorana fermions 

To study the emergence of Majorana fermions in the 
simplest possible setting, we now focus on the V — \ 
case and consider values of the chemical potential satis- 
fying < i.e. intersecting only the / = branch 
of the spectrum. The Hamiltonian for this branch then 
becomes — ks2 — ^J■ + ms^, where we have explicitly 
included the chemical potential term. The Bogoliubov- 
de Gennes Hamiltonian describing the proximity-induced 
superconducting order in the nanowire can be written, 
in the second-quantized notation, as H = X^fc^l^fc^fc 
with *fc = {fk,gkjik'9-kV and 

For the surface state, rj — 1/2 represents a T- invariant 
point at which /ilj, — hk- Therefore, Hk can be written 
as 

* - {X -t) ■ 

In the following, we consider the simplest s-wave pairing 
potential = Ao«S2, with Aq a (complex) constant 
order parameter, which corresponds to the pairing term 

^o{fk9-k ~ 9kf^k)- It useful to note that this form of 
Afc actually implies a vortex in the SC order parameter, 
as can be seen by transforming Jik back into the original 
electron basis, i.e. undoing the transformation indicated 
in Eq. ^ . The phase of the order parameter in this basis 
winds by 2t: on going around the cylinder as required in 
the presence of the applied magnetic field whose total 
fiux is $o/2. 

Introducing Pauli matrices Tq in the Nambu space and 
assuming Aq real, we can write 

Hk ^ T3{ks2 - n) + T^ssm - T2S2A0. (9) 

(Here we have again taken v = h — 1.) The spectrum for 
this Hamiltonian is Ek = ±{k'^ + fi'^ + m'^ + Al±2{k'^fj,'^ + 
^^w? + to^Aq)^/^)^/^. We now consider a special case 
when /X = 0. The Hamiltonian simplifies, T-Lk = T3S2k + 
TsSsiTi — T2S2A.0 and the spectrum assumes a simple and 
suggestive form 

Ek = ±V^-2 + (m±Ao)2. (10) 

We observe that the spectrum is fully gapped in the pres- 
ence of either SC or magnetic order but has a gapless 
branch when m = ±Ao. Thus, we expect a topological 
phase transition at this point. Consequently, we expect 
gapless modes to exist at an interface between SC and 
magnetic domains in a wire. 
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FIG. 4: A convenient possible choice for the SC/magnetic 
domain wall at 2 = 0. A{z) is the SC order parameter and 
m(z) the magnetic order parameter. A dashed line shows the 
zero-mode solution jj,{z) for this domain wall. This particular 
choice of boundary conditions can be used to show Majorana 
fermions occur at the ends of the TI nanowire irrespective of 
precise boundary conditions. 



Now consider spatially varying m{z) and A(z) such 
that to(0) = A(0) as sketched in Fig. [i) With these 
choices for the order parameters, we expect the spectrum 
to be gapped far away from the domain wall, but we 
expect gapless modes localized near z = 0. To determine 
if there are any fermionic zero modes, we rotate Hk in 
s — r space so that the rotated Hamiltonian is completely 
off-diagonal. That is, we work with Hk = UHkU~^, 
where 

U = e-'T^^e-'^^^ (11) 
Then Hk = TiS2k + Tisim — T2S2A., so Hk is of the form 

where Dk = + sim + is2A. 

We now replace k — > —id^ and look for solutions ^{z) 
satisfying 



(13) 



Taking ^'(z) = ("01 j "02 7 "03 7 "04) and reinstating v, Eq. 



( 13 1 yields four independent equations: 



{+vd^ + m + A)^pl = 
{-vd, + m- A)V'2 = 
(-l-i;^^ + m - A)V'3 = 
{-vd;, + m + A)ipi = 



(14) 
(15) 
(16) 
(17) 



Here we have suppressed the z dependence. The solution 
u{z) of an equation of the form 

[vd,+uj{z)]u,^0, (18) 

a Jackiw-Rossi zero mod^^, can be written as 

u{z)=uoe^-fod-'-'i-')_ (19) 



This solution is normalizable provided that (jj{z) has a 
soliton profile, i.e. is proportional to sgn(z) for large \z\. 
According to our assumptions, m(z) + A(z) > for all 
values of z, so there is no normalizable solution for -01 or 
-04. With V > 0, there is also no normalizable solution 
for ^02, but there is one for ip^: 



03(z) = Woe" 



■ dz'{m{z')^A{z')) 



(20) 



(For u < 0, "02 would be the normalizable solution in- 
stead.) Thus, for A{z), m{z) as given in Fig. |4j our 
Hamiltonian has a single zero-mode solution of the form 
'I'o = (0,0, l,0)'^u(z) localized near the domain wall at 
z — 0. This solution is valid as long as w > and 
■m{z) — A(z) — ^ ibconst for z — >• ±00. To see if the 
zero mode ^o(^) corresponds to a Majorana fermion, we 
undo the unitary rotation and inspect the corresponding 
solution ^o(z) = C^^^^o(2), which is 



■^^{z) = \[l -1,1, -If u{z). 



(21) 



In second quantization, the field operator destroying the 
particle in the state '^o{z) is 

V-o = i y dzu{z) [f{z) - g{z) + /t(z) - g^z)] . (22) 

where f{z), g{z) are real-space versions of the fk, gk 
operators in ^k ■ Since u{z) is real, it holds that 0^ — ipo, 
so represents a Majorana fermion. 

With a few additional observations, the above calcu- 
lation can be used to show that an additional unpaired 
Majorana mode exists at the SC end of the wire irre- 
spective of boundary conditions. First, recall that in a 
finite system Majoranas always come in pairs, since they 
are formed from ordinary fermion^. This second Ma- 
jorana fermion, being a zero-mode, cannot not exist in 
the nanowire bulk where the spectrum is gapped. It can- 
not exist at the magnetic end because the magnetic or- 
der does not support the requisite particle-hole mixing. 
The second MF must therefore be at the SC end of the 
nanowire, irrespective of the exact boundary condition. 
From this, we can argue that the special conditions used 
to establish the existence of unpaired MFs in the device 
are unnecessary: The zero-modes in fact exist in the de- 
vice under generic boundary conditions as confirmed by 
explicit numerical study using a lattice model, discussed 
in section HI. A specific example of Majorana end-states 
obtained in such a lattice calculation under general con- 
ditions is given in Sec. IV. A below. 



C. Energy gap protecting Majorana zero-modes 

As mentioned irP^, in order to detect and manipulate 
MFs under experimentally accessible conditions it is cru- 
cial that they are protected from all other excitations by a 
gap. The latter is often refered to as a 'minigap' because 
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typically there will be other excitations inside the bulk 
gap. We study the minigap in this TI nanowire-based 
device both analytically and numerically. 

In this section we estimate the minigap for the super- 
conducting TI nanowire using the analytical low-energy 
theory. Specifically, we wish to fin d th e lowest non- 
zero eigenvalue of T-Lk defined in Eq. (12 1. We start by 
squaring the Hamiltonian. We find, with k — >■ —idz and 



DlDk 



D^D, 



D^D = dl + [A'(z) - s^Tn!{z)\ + [A(z) - s^m{z)f . (23) 

The two independent equations for S3 = ±1 in D can 
more conveniently be written as {D'^U)j^ and (D^D)^, 
where 

{D^D)^^dl + [A'{z)Tm'{z)] + [A{z)Tm{z)f ■ (24) 

To find the energy of the first excited state, we look for 
solutions ip satisfying 



(25) 



where h = 'H? and e > 0. We consider m{z), A(z) such 
that A(z) -I- m(z) =:const for each z, and A(z) — m(z) — 
f{z) having a soliton profile, e.g. we may take f{z) = 
Aotanh(2/^), as shown in Fig. Wl Then {D^D)_ yields 
no bound states. (D'''D)+, however, has the form, with 
velocity v restored. 



v'dl+vf'{z) + f{z) 



(26) 



For bound-state energies much less than Aq, f{z) can be 
approximated as linear in the vicinity of = 0. With 
f{z) ~ — Ao|, where ^ is the length scale over which the 
SC order parameter varies near the domain wall, we then 
have 



2a2 



-v'd 



vAo 



Ao 



(27) 



This is the Hamiltonian for the harmonic oscillator with 

... 2 

the identification 



2 



if)' 



and h'^uj'^ 



4w^^. Therefore, allowed eigenenergies of {D'^D). 
bounded above by A§ are 

e„ = hw{n -I- - 



uAo 2hvAQ 



(28) 



where n is any non-negative integer. The energy spec- 
trum of T-L in this approximation is then 



2hvAQ 



(29) 



Since ^ is the length scale over which the SC order pa- 
rameter varies near the wire end, it is at most the SC 
coherence length The minimum energy of the first 

excited state Ei is then 



El = Aov^, 



(30) 



which is already greater than Aq. Therefore, there are 
no excited states where the linear approximation holds. 
There can be some at energies close to Aq and this is 
consistent with numerical results presented irP^. We thus 
conclude that the minigap amplitude is close to Aq in this 
case. 

We also note that the calculation presented above is 
valid in the special case /i = 0. For non-zero chemical 
potential the situation is more complicated and we are 
not able to find a simple analytic solution for the excited 
states in this case. Since the density of states of the 
underlying Dirac semimetal grows with increasing energy 
we expect there to be more low-lying excited states when 
/i 7^ and thus reduced minigap. This expectation is 
indeed confirmed by our numerical simulations discussed 
below. 



D. Majorana state in a finite-wire configuration 

So far in our analytical calculation we have shown that 
a Majorana state exists in a TI wire at the interface be- 
tween SC and magnetic domains. A more realistic situ- 
ation from an experimental point of view is to consider 
a finite-length wire located on top of a s-wave supercon- 
ductor. One would expect the Majorana zero modes to 
live close to the ends of the wire where the SC order 
parameter vanishes and the magnetic flux pierces of the 
surface of the wire. Our numerical results below confirm 
this intuition but here we briefly show that one can also 
analytically prove this for a simple wire configuration and 
find the Majorana state. 

Consider one end of a finite wire with a configuration 
as shown in Fig. ([s]): a cylinder with the sharp edge 
smoothed out. For convenience we think of the result- 
ing surface separated into three regions. As before we 
consider a uniform magnetic flux with ry = 0.5 through 
the wire but we assume that it only pierces region C of 
the surface. The form of the Hamiltonian is different in 
each region as one would expect due to the curvature ef- 
fects and is discussed in detail in the Appendix B. The 
regions have been chosen in such a way that the first spa- 
tial derivative operator that appears in the Dirac Hamil- 
tonian is continuous everywhere on the wire including 
the boundaries between the regions. This way it is legiti- 
mate to use Dirac hamiltonian Eq.Q to study the surface 
states since the normal unit vector to the surface of the 
TI is well-defined everywhere on the wire. Note that this 
would not be the case if we considered sharp edges. We 
can exploit the azimuthal symmetry of this configuration 
and perform the unitary transformation defined by Eq. 
Q. The full BdG Hamiltonian for the I = branch then 
reads 



H 



(31) 



where 'H(C) is a 4 x 4 matrix given hy {h = v ~ 1) 
m) = [^S2dc + giOsi + m(C)s3 + A(C)] - A(C)t2S2 
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FIG. 6: The probability density of the Majorana state (up 
to a normahzation constant) close to one end of the wire for 
Ao — hv/Ro and ao = 0.17?o- Note that l^l'^ peaks at the 
boundary between toroidal region and the disk-like end. It 
decays exponentially into the bulk of the wire (region A as 
shown in Fig. ([5|) and decays with a power law behaviour 
toward the centre of the disk (region C). 



FIG. 5: The schematics of a finite-length wire on a supercon- 
ducting substrate. The surface of the wire has been divided 
into three regions A, B and C (bottom panel) according to 
the form of the Hamiltonian. The top panel details the as- 
sumed shape of the wire end with various quantities used in 
the text indicated. 



The Hamiltonian is a function of the length Q which 
parametrizes the geodesic curves that connect two end 
points on each section of the surface. We choose C, to 
be zero at the end of the cyhnder (i.e. at the boundary 
between regions A and B as shown in Fig. Note 
that on the cylinder it is equivalent to the z variable we 
used before. In general we have 



z (z < 0) 



CeA 



C,= {aoe (0 61 7r/2) C e ^ (32) 

.i?o + ao^/2-p (Os^p<i?o) CeC 

where Oq , Rq are the radii of the connecting torus and the 
cyhnder (ag <C Ro) respectively. g[C,) and \{C) are two 
functions that arise due to the curvature effects. They 
are defined as 



5(C) 



1 







1 - [p/Rof C e C 



and 



1 (0 CeAC 
A(C) = X 



(33) 



(34) 



Assuming that the magnetic flux is narrower than the 
wire cylindrical shaft we can neglect the Zeeman field in 



all regions except region C . Therefore we consider the 
following profile for it 



m(C) 




(35) 



As mentioned previously, we expect MF to exist near the 
end of the wire irrespective of the details of the boundary 
condition. For simplicity, therefore, A((^) is assumed to 
be nonzero and uniform only in the range of C, parametriz- 
ing the cylindrical shaft of the wire and the rounded edge 
(C e A,B). In region C we assume A((^) = in accord 
with the intuition that the SC order will be suppressed 
here due to the magnetic field piercing the surface and 
the presence of the vortex. Finally, we consider a very 
long wire so that we can assume that the overlap between 
localized states at the two ends is negligible. This way 
we can look for solutions with zero energy and safely ex- 
clude the solutions that grow exponentially towards the 
other end. Thus, we seek the solutions of the following 
equation 



^(C)*(C) = 



(36) 



and investigate whether there is a solution (localized near 



the end) which satisfies the Majorana condition = 'I'o 



The strategy for solving Eq. (36) is standard: we find 



the general solutions in the three regions A, B and C, we 
match their wavefunctions at the boundaries and finally 
select the physical (normalizable) solution. The technical 
details are given in the Appendix B. Up to a normaliza- 
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tion constant the solution is given by 
where the wave functions are defined as 



(37) 



(38) 



COS : 



COS : 



sm 



ysin 



(39) 



I -cosy 



Assuming that the Zeeman term is negligible and con- 
sidering only the solution which is well defined (i.e. no 
ambiguity in the phase at the apex which implies that 
the wave function should vanish at p = 0) we obtain 



*c(p) 



2i?f 



exp 



1 + 2Aoao7r) 



(40) 



The probability density of the Majorana state is shown 
in Fig. ([5]). Note that although the components of the 
wavefunction change in the B region the absolute value 
follows the same behaviour as a function of length (. The 
peak of the wavefunction is at the boundary between re- 
gion B and C. This way we obtain a solution to the 
Bogoliubov-de Gennes Hamiltonian which is real (up to 
an overall phase) and it is associated with a zero energy 
eigenvalue for a semi-infinite wire. The wavefunction is 
localized at the end of the wire in agreement with the nu- 
merical simulation and satisfies the Majorana condition 



IV. RESULTS ON STABILITY OF MAJORANA 
FERMIONS 

A. Lattice model 



with Mk = e — 2t '^^^ ^a- Here ctq, represent the Pauli 
matrices acting in the space of two independent orbitals 
per lattice site. For A, > and 2t < e < 6t the sys- 
tem described by Hamiltonian (41) is a TI in Z2 class 
(1;000), i.e. a strong topological insulator. The magnetic 
field enters through the Peierls substitution, replacing all 
hopping amplitudes as — >■ tij cxp [— (27rj/$o) // A • dl] 
and the Zeeman term — g^^B • s/2 where /is = eH/2meC 
is the Bohr magneton. In the SC state the BdG Hamil- 
tonian takes the form of Eq. ([t]) with Ak = Ao«S2 de- 
scribing on-site spin singlet pairing. 

In the subsequent calculations we consider the above 
Hamiltonian on the real-space cubic lattice and in various 
wire geometries with rectangular cross sections and both 
periodic and open boundary conditions along the length 
of the wire. We find eigenstates and energy eigenvalues 
by the exact numerical diagonalization using standard 
LAPACK routines and by sparse matrix techniques in 
cases where only low- lying states are of interest. Unless 
explicitly stated otherwise we use the following set of 
model parameters in our subsequent calculations: A^ — 
2A, t = X, e = 4A. This places our model into Z2 class 
(1;000) and with A = ISOmeV produces a bulk bandgap 
of 300meV, as in Bi2Se3 crystals. 

As an example of a calculation based on the lattice 
model we show in Fig. [7] the probablity density of Majo- 
rana end-states in a wire 100 lattice spacings long. We 
note that in any finite-length wire there will always be 
an exponentially small overlap between the two Majo- 
rana end-states. Such an overlap leads to the hybridiza- 
tion and a small non-zero energy SE for the combined 
fermionic state which shows probablity density equally 
split between the two ends of the wire. An (unphysi- 
cal) state with equal probablity density exists at energy 
—SE. Fig. [7] shows how the Majorana end-states can be 
constructed by taking the appropriate linear superposi- 
tions of the above eigenstates. It is to be noted that for 
a finite-length wire the Majorana end-states are not true 
eigenstates of the system; they will mix on a time scale 




We establish the stability of Majorana fermions in the 
nanowire through a combination of additional analytical 
insights and numerical studies using the same concrete 
lattice model for the Bi2Se3 family of materials^ 
given by Fu and Berg^* regularized on a simple cubic 
lattice. This model is defined by a /c-space Hamiltonian 

/ik = AfkCTi + Aa-3(s2 ainkx-si smky) + \zcr2 sinfc^, (41) 



FIG. 7: Probability densities of Majorana end-states. Pa/b{z) 
represents the particle component of the wavefunction associ- 
ated with zkEo summed over x-y coordinates. Pa±b represent 
the even/odd superpositions of these wavefunctions. The wire 
is 6 lattice sites wide in both the x and y directions and 100 
lattice sites long in the z direction. 77 — 0.49 in units of 
the fundamental flux quantum, and the chemical potential 
H = 0.09A, where A = 150 meV. 
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h/5E which is, however, exponentially long in the wire 
length L. 



B. The Majorana number and the topological 
phase diagram 

The existence of MFs at the ends of the wire depends 
on whether or not the wire is in the topological phase. 
For a ID system the presence of the topological phase 
is indicated by Kitaev's Majorana mimhe^^ A4{H). In 
Ref."^ we have computed M{H) in the limit A — 
where it reflects simply the parity of the number of the 
Fermi points of the underlying normal state in the right 
half of the Brillouin zone. The resulting phase diagram 
in the rj-fj, plane, for /i inside the bulk bandgap, consists 
of diamond-shaped topological regions indicated in Fig. 
[sjja). In the limit A — > the individual diamonds just 
touch at their apices yielding a continuous topological 
phase for a specific value of magnetic flux rj close to | 
and for all values of fj, inside the band gap. This feature 
underlies the key advantage of the present setup: the 
chemical potential does not require fine tuning. How- 
ever, it is also true that when passing between individ- 
ual diamond-shaped regions, the system gets arbitrarilly 
close to the phase boundary and one thus expect topo- 
logical order to be rather fragile in these regions. On the 
other hand, away from the A — >■ limit one intuitively 
expects the topological phase to become more robust and 
indeed this was suggested by the numerical results pre- 
sented in Ref.E^. In the following we shall elucidate this 
point and show that indeed for Aq > the topological 
phase becomes a compact region in the rj-fj, plane. 

We consider a general definition of the Majorana 
numbeil^^ 



M{H) 



sgn 



F{B{k = 0) 



sgn 



a 



(42) 



where B{k) is the position-space Hamiltonian of the 
infinite-length, lattice-model TI wire written in terms of 
Majorana fermion operators and Fourier-transformed in 
the ^-direction, a is the lattice-spacing, implying that 
B{k) is evaluated here at k = and at the edge of the 
Brillouin zone. Eq. ( 42 ) simplifies to the previously- given 
definition irP^when^Aj is sufficiently small. 

We consider the phase diagram in two limiting situa- 
tions: When A winds in phase counterclockwise by 2tt 
around the circumference of the wire, corresponding to 
a vortex present along the length of the TI nanowire, 
and when A does not wind in phase, meaning there is no 
vortex in the system. These two situations are expected 
to represent the ground state of the system for the mag- 
netic flux close to $0/2 and 0, respectively. The precise 
value of the flux at which the vortex enters will depend 
on details but we show below that, remarkably, the topo- 
logical phase diagram is fairly insensitive to the presence 
or absence of the vortex. 

Fig. [Sjdisplays the phase diagram of a rectangular wire 
with a 10 X 10 cross section based on Eq. (42 1. As |A| 



is increased from zero, we see that the boundary of the 
region corresponding to the topological phase smoothes 
out, with the region centered close to 77 = 0.5 and extend- 
ing through all values of the chemical potential inside 
the bulk band gap and also up into the bulk conduction 
band. We remark that the topological phase here is cen- 
tered near a value of the magnetic flux that somewhat 
exceeds $o/2. This is because the surface state pene- 
trates slightly into the bulk of the wire and thus encloses 
somewhat smaller amount of flux than the geometric sur- 
face area of the wire. For thicker wires this shift will be 
negligible. 

It is also interesting to note that according to Fig. [8] 
the topological phase persists for /i well into the conduc- 
tion band (as well as the valence band, which is not ex- 
plicitly shown). This finding is potentially important in 
view of the fact that most TI crystals naturally grow with 
the chemical potential pinned inside the bulk conduction 
or valence band. We will show below, however, that al- 
though MFs indeed appear in this regime, the minigap 
protecting them quickly collapses as fi moves deeper into 
the bulk band. Heuristically, one can understand this as 
follows. With /i inside the conduction band the bulk of 
the nanowire becomes metallic. In the presence of SC 
order and with the magnetic field applied along its axis 
there will be a vortex line running along its center. Such 
a vortex line will host low-energy core states with the 
characteristic energy scale A'^ /Ep which quickly becomes 
very small as Ep increases (here Ep is the Fermi energy 
measured relative to the bottom of the bulk conduction 
band). By contrast when the chemical potential lies in- 
side the bandgap the bulk of the wire remains insulating 
and does not contribute any low-energy states. 

From study of the phase diagram Fig. [8] at A > 0, 
we begin to understand the robustness of the Majorana 
bound states. Consider first the effect of non-magnetic 
disorder, modeled by introducing a spatially fluctuating 
component of the chemical potential fj, ^ fi + Sfxir). As- 
sume also that 5fi{r) is slowly varying in space so that 
only variation along the z-direction are meaningful and 
fi{z) can be thought of as deflning a phase of the wire 
in the vicinity of the coordinate z. If the average chem- 
ical potential p, and the flux 77 are such that the system 
starts deep inside the topological phase then it is clear 
that fluctuations in dfi will not drive the system out of the 
topological phase unless they exceed the bulk bandgap. 
Similarly, fluctuations in the total magnetic flux 77, which 
can occur e.g. in a wire with a non-uniform cross section, 
will not drive the system out of the topological phase un- 
less they reach a significant fraction of $0/2. We demon- 
strate below by explicit inclusion of disorder in the lattice 
model that the heuristic argument given above remains 
valid even when disorder potential varies rapidly on the 
lattice scale. 

The smoothing out of the topological region's bound- 
ary as |A| is increased can be understood by studying 
the low-energy analytical theory again. We start from 
the Hamiltonian in Eq. (IS]) and let n = m = 0, study- 
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FIG. 8: Phase diagrams computed from generalized Majorana number for the infinitely-long, clean wire with a 10 lattice site 
by 10 lattice site cross-section, (a), (b), and (c) are each for a system with a vortex, while the system with phase diagram (d) 
lacks a vortex. |A| is set to and 0.04A in (a), (b), respectively, and 0.08A in both (c) and (d). Chemical potential fi is in 
units of A = 150 meV and rj is in units of the fundamental flux quantum. Blue and pink regions were created by computing 
the generalized Majorana number A4 in steps of at most A-q — 0.02 and — 0.02, colouring these points blue {A4 = 1, 
non-topological phase, no Majorana zero- modes in system) or pink {M = — 1, topological phase, Majorana zero-modes present 
in system), respectively, and enlarging these data points to create regions of solid colour. As well, the phase boundaries were 
computed (white lines) via a more efficient algorithm, with error bars at most the size of the symbols. White lines extend only 
up to /X = 1 as at larger fi the topological phase regions break up into small domains and the algorithm used to compute the 
phase boundary is only effective for large, simply-connected regions of the phase diagram. 



ing how the phase diagram changes for this value of the 
chemical potential as |A| is increased from zero. If we 
now assume that rj deviates from 1/2 by a small amount, 
i.e. rj = 1/2 + 6'!] then the spectrum for the I — branch 
can be written as 



± 



± Ao 



1/2 



(43) 



We know that the system will be in the topological 
phase when Si] = and /i = 0*. To understand the 
smoothing out of the topological region's boundary, we 
identify when the spectral gap closes for non-zero Aq 
as a function of Srj, since closing of the gap signifies a 
phase transition. We notice that the gap in Eq. ( 43 ) re- 
mains as Sr] is moved away from until Si] — zLRAq/vH, 
where we have restored proper units. Therefore, we see 
that at = 0, for nonzero Aq, the topological phase 
has widened from a point at 77 = 1/2 to an inter- 
val {1/2 - RAo/vh,l/2 + RAo/vh), as observed in the 
phase diagrams computed numerically using the more 



general definition of the Majorana number. 

The absence or presence of a vortex in the TI nanowire 
makes negligible differences to the phase diagrams as seen 
by comparing Figs.jsjjc) andjsjjd). However, the presence 
or absence of a vortex does have a pronounced effect on 
the quasiparticle excitation gap Aoxc (shown in Fig. [9]) 
in the infinitely-long wire. From Fig. [9] we see that 
Aexc remains close in magnitude to | A| up until fi reaches 
the bottom of the bulk conduction band if the SC order 
parameter winds counter-clockwise in phase by 27r, while 
without a vortex Acxc can be seen to quickly decay with 
increasing /x. Clearly, near 77 = 1/2 the former represents 
a more physical situation. 

For experimental realization of the device, it is im- 
portant to know how large the minigap - the energy of 
the lowest non-zero mode coming from bound states at 
the ends of the wire which are absent when the wire is 
infinitely long - is as well as studying the quasiparticle 
excitation gap, which is the energy of the lowest non- 
zero mode for the infinitely-long wire. If we study this 
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FIG. 9: Quasiparticle excitation gap Aexc of an infinitely-long 
wire without disorder with a 14 lattice site by 14 lattice site 
cross-section as a function of chemical potential jj, with vor- 
tex present (black circles) and vortex absent (red diamonds). 
Aexc and fi are expressed in units of A = 150 meV. Here 
|A| = 0.08A. 

minigap as a function of chemical potential with a vor- 
tex present, shown in Fig. [TOj we see that, for fi close to 
zero, the minigap starts with an amplitude close to the 
SC gap A, in accord with the analytical theory presented 
in Sec. III.B. The minigap then continuously decreases 
with increasing /i, retaining a respectable value ~ O.IA 
at the edge of the bulk conduction band at = 1. As 
mentioned previously, the minigap then quickly collapses 
as the bulk bands are populated but nevertheless persists 
over a non-zero range of fi > 1. 

C. Robustness of Majorana fermions against 
disorder 

As mentioned previously, we expect robustness of the 
Majorana end states with respect to non-magnetic disor- 
der in the proposed device. There are two related but log- 
ically separate issues that pertain to this problem. First, 
we must ensure that SC order induced in the wire is itself 
robust against non-magnetic disorder. Since the device 
can be operated at (or very close to) the time-reversal in- 
variant point and since we consider a spin-singlet s-wave 
SC order we expect this to be the case on the basis of 
Anderson's theorem. Below, we illustrate this aspect of 
the robustness by performing a self-consistent numerical 
calculation on our model in the presence of disorder. Sec- 
ond, we must show that Majorana end-states themselves 
remain stable in the presence of disorder. To some extent 
this already follows from our arguments in the previous 
subsection based on the study of the topological phase di- 
agram. Also, stability of MFs follows from the stability of 
the SC phase in the bulk of the wire argued above. Never- 



theless, to address this question directly, we perform ex- 
plicit numerical calculations for open-ended wires in the 
presence of non-magnetic disorder and in various phys- 
ical regimes. These calculations confirm the expected 
robustness of MFs and yield additional insights into the 
quantitave aspects of this robustness; specifically they 
provide information about the minigap magnitude and 
the mechanism by which MFs are eventually destroyed 
in the strong-disorder limit. 

To study these questions, we add a term corresponding 
to disorder in the on-site potential to the lattice Hamil- 
tonian, Hq, of the clean system. The Hamiltonian for the 
system with disorder, -ffdis, is therefore 

Hdis = Ho+Y^ U,c\^c,^> , (44) 

ia 

where C/^, the random potential at lattice site i, is as- 
signed a value from a uniform, random distribution rang- 
ing from ^ to ^. 

As a first step we compute the magnitude of the SC 
order parameter self-consistently as described irP^l for 
different disorder strengths and with periodic boundary 
conditions along z, i.e. no Majorana end-states. This 
calculation assumes the existence of an intrinsic pairing 
potential V in the wire and is therefore, strictly speak- 
ing, not directly relevant to the proximity-induced SC 
state discussed in the rest of this paper. It neverthe- 
less illustrates very nicely the robustness of the SC order 
with respect to disorder. The key point is that one would 
expect SC order to be even more robust when induced 
by proximity to the bulk superconductor. Fig. |11| shows 
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FIG. 10: Three lowest energy eigenvalues E2, E\, and Eo of 
the finite-length wire, without disorder, with 14 by 14 by 100 
lattice sites as a function of chemical potential /x with vortex 
present with |A| = 0.08A. l^exc and ^ are expressed in units 
of A = 150 meV. Eigenvalues were computed via Lanczos 
method and failed to converge for ^ > 1.05, resulting in a 
non-physical spike to the non-convergent next data point. 
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lowering of the excited states at some critical value of 
disorder strength Uc- It is interesting to note that Uc is 
rather large, being expressed in units of A = 150 meV, 
exceeding the TI bulk bandgap by more than a factor of 
5. Furthermore, the MFs are robust against disorder at a 
wide range of average chemical potential in the bulk band 
gap, although the minigap is largest for values close to the 
middle of the bulk band gap and smallest for values near 
the edge of the bulk band gap. We note that the minigap 
is roughly the same for mean chemical potential values 
of 0.09A, 0.4A, and 0.8A if ^ = A, suggesting disorder 
stabilizes the minigap in this regime. It is also interesting 
to note that for two larger values of fj, disorder initially 
increases the minigap thus making the topological phase 
more robust. 



FIG. 11: Mean superconducting order parameter magni- 
tude |A|avg for three different values of the pairing potential 
{V = 1.310, V = 1.295, and V = 1.290) as a function of dis- 



order strength 



is averaged over every lattice site in 



a 6x6x6 lattice site nanowire with periodic boundary condi- 
tions in the £-direction with random disorder in the chemical 
potential and also further averaged over 10 such randomly- 
disordered nanowires. The mean chemical potential for all 
data points is /Xavg = 0.09A. The self-consistent calculation 
for each disordered nanowire proceeded until the maximum 
difference in the superconducting order parameter magnitude 
between the final iteration and the next-to-Iast iteration of 
the calculation at any lattice site was 0.001 A. 



mean SC order parameter magnitude vs. ^ computed for 
different values of the pairing potential V. We observe 
that |A|avg first decreases slightly with increasing U, but 
at larger U the mean SC order parameter increases in 
magnitude. We attribute this increase to the buildup of 
the normal density of states at the Fermi level, N{fi), 
due to disorder. Such a buildup is known to occur in 
other 2D systems with a Dirac spectrunPSl and it should 
increase |A| according to the standard BCS formulaPSl 



-i/7V(p)y 



(45) 



Here ujc and V are constants. These results suggest the 
SC order parameter is not only quite robust against non- 
magnetic disorder but the latter can acually enhance it 
when the chemical potential is close to the Dirac point. 

To address the robustness of Majorana end-states we 
now study the finite-length wires. Using sparse matrix 
techniques, we solved for the average values of the three 
lowest, positive eigenvalues of the spectrum, and plotted 
these for a range of U, as shown in Fig. [T2j These 
calculations are performed for a constant value of the 
SC gap, having previously established that disorder has 
only a weak effect on the latter. We see that the energy 
of the Majorana bound state remains very close to zero, 
with no observable fluctuations. The topological SC is 
eventually destroyed by the collapse of the minigap, i.e. 



V. CONCLUSION AND DISCUSSION 

The main goal of this work was to study the stabil- 
ity of Majorana zero-modes in the proposal introduced 



irP^, which consists of a TI nanowire, proximity-coupled 
to a bulk s-wave superconductor, with a weak magnetic 
field applied along the nanowire's length. After reviewing 
the literature to illustrate the importance of identifying 
a device in which MFs emerge under a wide range of ac- 
cessible conditions, and reviewing the theory behind the 
proposed device in greater detail than possible u^^, we 
studied the robustness of the topological superconductor 
phase of the device against disorder. 

As a first step we computed the topological phase di- 
agram of the TI nanowire numerically in a semi-realistic 
lattice model. Using a general definition of Kitaev's Ma- 
jorana number Ai{H) we were able to show that for 
non-zero values of SC order parameter A the topologi- 
cal phase forms a set of compact columnar regions in the 
plane spanned by the magnetic fiux rj = ^/^o smd the 
chemical potential /i, centered around half-integer values 
of v and covering about 50% of the phase diagram (see 
Fig. 8]). This form of the phase diagram confers two prin- 
cipal advantages of our proposed device as regards future 
experimental realizations and potential applications: (i) 
unlike in the semiconductor wire realization^! where sig- 
nificant fine-tuning of the chemical potential is required 
to reach the topological phase, our proposed device pro- 
duces Majorana end-states for /i anywhere inside the bulk 
bandgap of the TI; and (ii) if the average chemical po- 
tential of the nanowire is in the bulk band gap, we can 
expect the entire nanowire to remain in the topological 
phase even for large disorder strengths, as even then any 
local chemical potential value will remain in the topo- 
logical region of the phase diagram. We also observe 
from the phase diagram Fig. |8] that the topological phase 
persists over a wide range in magnetic flux through the 
cross-section of the nanowire. This feature is less critical 
because the magnetic field can be easily tuned in a lab- 
oratory. Nevertheless understanding the magnetic phase 
boundary is very useful since changing the magnetic field 
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FIG. 12: Three lowest positive energy eigenvalues Eq, Ei, and 
E2 of a finite-length TI nanowire model with 6 by 6 by 100 
lattice sites obtained by the Lanczos method as a function 
of disorder strength ^ for mean chemical potential values of 
H = 0.09 (a), fi = 0.4 (b), and /x = 0.8A. |Aj = 0.08A for each. 
Y and energy E are expressed in units of A = 150 meV. The 
error bars reflect averaging over 10 independent realizations 
of the random potential. 

strength can be used to tune the wire between topological 
and normal phases. 

To explicitly ascertain the robustness of the MF with 
respect to thermal fluctuations and non-magnetic disor- 
der we studied the system's minigap by both analytic 
and numerical techniques. Minigap, defined as the small- 
est non-zero energy level in the system, can be taken as 



a good proxy for the robustness of the quantum infor- 
mation encoded in MFs since uncontrolled excitations 
of quasiparticles out of the ground state into the low- 
lying excited states would obviously spoil such encoding. 
Also, large values of the minigap can aid experimental 
detection of the Majorana zero mode using various spec- 
troscopic techniques. We use the low-energy, analytical 
theory oP^ to first show that excited states should be 
close in energy to the magnitude of the superconducting 
gap when the chemical potential is close to the middle 
of the bulk band gap of the TI nanowire, indicating Ma- 
jorana zero-modes in the device should be protected by 
a sizeable minigap. We then employ the lattice model 
and compute the three lowest eigenenergies of the TI 
nanowire with disorder to show that the minigap remains 
significant even for the disorder strength considerably ex- 
ceeding the bulk band gap of the TI. Interestingly, we 
find that disorder strength comparable to the magnitude 
of the bulk bandgap also appears to stabilize the minigap 
at a sizeable value over changes in the average chemical 
potential in the nanowire, which might be useful in future 
applications of the device. 

Stability of MFs in our proposed device also follows 
from the general periodic classification of topological in- 
sulators and superconductors given by Schnyder et alW^ 
and by Kitae\i22i. According to this classification Majo- 
rana bound states may appear in one-dimensional sys- 
tems of the D and Dili symmetry classes. The former 
corresponds to the superconducting state with broken 
time-reversal and spin symmetries while the latter is a 
superconductor with only spin symmetry broken. We 
found that MFs are most stable in our device with ex- 
actly half fiux quantum, corresponding to T-preserving 
Dili symmetry class. In this class all the states are dou- 
bly degenerate according to the Kramers theorem. It is 
thus slightly counterintuitive that we obtain isolated non- 
degenerate MFs in this situation. To clarify this point 
we note that the device will be in the Dili class only 
when the wire is considered infinitely long (or else peri- 
odic along the z-direction) . MFs on the other hand occur 
only in a finite wire where T is explicitly broken at the 
wire ends by the magnetic fiux lines entering and exiting 
the wire. 

The above considerations also suggest stability of the 
MFs with respect to magnetic disorder which was not 
explicitly studied in this paper. For one, magnetic im- 
purities are much less dangerous for proximity-induced 
SC than intrinsic SC. Furthermore, magnetic impurities 
will not further violate the class D symmetry possessed 
by this system when T is broken. 

Before concluding we address the following question: 
Do MFs predicted to exist in our device obey non- 
Abelian exchage statistics? This property is obviously 
of paramount importance for any future application in 
quantum information processing. Alicea et al.'^^l clari- 
fied the sense in which MFs in ID wire networks exhibit 
non-Abelian statistics upon exchange, considering semi- 
conductor wires^'^'^^', and showed how particle exchanges 
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exchange statistics by virtue of being adiabatically con- 
nected to the Fu-Kane vortex^^ and, as illustrated above, 
can be easily manipulated by changing the direction of 
the applied magnetic field. On this basis we expect the 
device architecture discussed in this work to be useful for 
future experimental study of MFs and their their poten- 
tial applications. 



FIG. 13: Exchange of MFs in a trijunction device, a) The 
field B is tuned so that nanowire 1 has flux close to $o/2 
through its cross section and is thus in the topological phase 
with MFs localized near its ends. The flux through wires 2 
and 3 is down by the factor cos(27r/6) = 0.5 and they are thus 
in the trivial phase, b) Rotating the direction of B by 30° 
the flux through wires 1 and 2 becomes cos(27r/12)($o/2) ~ 
0.866("l>o/2) and is thus sufficiently close to "I>o/2 for them 
both to be in the topological phase according to the phase 
diagram in Fig. [S] As a result the MP previously located at 
the junction (red circle) has now moved to the other end of 
wire 2 as indicated. Continuing this process by rotating B 
further in 30° increments it is easy to map out the motion of 
MFs and conclude that after 180° rotation the system comes 
back to the original situation with MFs localized on wire 1 
but with their order exchanged as illustrated in panel c). 



can be effected in such a set ting. Although superficially 
similar to these modelJ^^'^ our proposal is more closely 
related to the original Fu-Kane vortex at the interface 
between a TI and an ordinary superconductor^. Indeed 
consider a thought experiment in which we take our wire 
and slowly increase its radius while simultaneously re- 
ducing the applied magnetic field so that the total flux 
through its cross section remains constant at $o/2. We 
also assume that all the surfaces of the resulting disk 
remain covered by a thin SC film. In the limit when 
the radius R becomes comparable to the wire length L 
we have a bulk disk-shaped TI covered with a SC film. 
The presence of $o/2 flux and the cylindrical symmetry 
dictates that a vortex must be present at the center of 
each of the flat disk surfaces. Such vortices will contain 
MFJ^ and will obey non-Abelian exchange statistics ac- 
cording to the standard arguments. MFs in our wires 
are thus adiabatically connected to those residing in the 
cores of Fu-Kane vortices and are therefore expected to 
obey the same non-Abelian exchange statistics when or- 
ganized into T-junctions or wire network geometrie^^. 
In Fig. [T3| we outline a simple protocol that implements 
an exchange of two MFs in a symmetric 'trijunction' de- 
vice formed by three nanowires joined at a single point. 
The exchange of two MFs, initially localized at the ends 
of one of the wires, is effected simply by a 180° rotation 
of the applied magnetic fleld in the plane of the device. 

We conclude that unpaired Majorana zero-modes are 
exceptionally stable in our proposed device, being present 
over a wide range in magnetic flux, chemical potential, 
and disorder strength, with disorder even being expected 
to stabilize the MFs to an extent. They obey non-Abelian 
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Appendix A: Surface state in a thin wire 

To study the surface state in a thin wire we employ 
a continuum version of the 3D lattice model defined in 



Eq. (41 1 . Focusing on the vicinity of the F-point we thus 



write the following 4x4 matrix Hamiltonian, 

hk = -icr3{s2TT.j, - SiTTy) + \z(T2k^ + <7im{r), (Al) 

where tTj = ^idj — Aj is the gauge-invariant momentum 
operator and m(r) is the gap at the F-point which we take 
to have a radial depence to model an interface between a 
TI (m > 0) and an ordinary insulator or vacuum (m < 0) . 
We also assume traslational invariance along the wire, 
thus keeping the quantum number (which we denote 
as k hereafter). 

We start by considering a 5- function flux A = ri^o{z x 
f)/27rr for which we can flnd an analytic solution. Trans- 
forming into the cylindrical coordinates we obtain 







e-^(_a, + 2 + ia^) 




We now perform a unitary transformation in the spin 
space defined by the matrix 



,e^'^' 



(A2) 



with / the integer angular momentum quantum number. 
This leads to a decoupling in the angular variable with 
the individual channels described by 



hki{r) 



132 



dr 



1 

2^ 



1 



1 + 



1 



+ (j2Xzk + aim{r). 



(A3) 



From our previous discussion of the surface state we 
expect the gapless mode to occur for Z = and rj = ^. 
Furthermore, the band crossing occurs at fc = so we 
focus on this value of k and study 



hoo{r) -10-382 \dr + + (7im{r) 



(A4) 
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By performing a 7r-rotation around ui and around si this 
Hamiltonian can be brought into an off-diagonal form, 



(A5) 

suitable for identifying the possible zero modes. We now 
describe the TI wire of radius R by the following config- 
uration of the gap function, 



m(r) = 



Too 



r < R, 
r > R, 



(A6) 



with Too, TOq positive constants. For a wire placed in 
vacuum we take TOq — oo. With this gap configuration 
it is easy to see that exactly two zero modes exist, 

^aW =XaC^/-e'"°^ r<R, a = 1,2, (A7) 

V r 



which we treat in the first order perturbation theory. It is 
straightforward to evaluate the requisite matrix element 
{tpi\Sh\ilj2)- For 7? = ^ and in the limit of C <C i? this 
leads to a correction to the energy of the form 



SE . ±iTOo 



(AlO) 



The energy splitting is seen to be proportional to the ra- 
tio between the cross-sectional areas of the surface wave- 
function and of the cylinder. Thus, we arrive at the con- 
clusion that the gap in the spectrum is indeed propor- 
tional to the amount of the T-breaking in the system as 
measured by the exposure of the surface wavefunction to 
the magnetic field. 



Appendix B: Dirac Hamiltonian for a finite wire 
configuration 



with C ~ \J 2mQ I {e^™-'-^^' — \)R the normalization con- 
stant and xi = (0,1,0,0)'^, X2 (0,0,1,0)'^. Thus, 
we conclude that exact zero modes exist for any radius 
R in the presence of a 5-function flux. This implies 
that the gapless surface modes also persist in this case 
even for a thin wire. We note however that only when 
i? ^ ^ = I/toq are the zero-mode wavefunctions local- 
ized near the surface; for a thin wire they permeate the 
entire bulk of the wire. 

The above result can be understood based on a sim- 
ple general argument. With the (5- function half-flux the 
3D Hamiltonian (All remains T-invariant. Therefore, 
Kramers theorem protects the band crossing at fc = 0, 
and the surface mode must remain gapless for arbitrary 
R. This understanding suggests that the degeneracy at 
fc = will be split by a more general flux distribution (e.g. 
that corresponding to a uniform B field), and the size 
of the gap will then refiect the 'amount' of T-breaking 
present in the system. Unfortunately, we were unable to 
find an exact solution for a more generic flux distribution. 

To test the above hypothesis we therefore pro ceed as 
follows. We start from the exact solution (A7| of the 
Hamiltonian /loo and add to it magnetic field (5B(r) as a 
perturbation. For simplicity we take the form 5B(r) = 
z(ar ~ b) with constants a and b chosen so that the total 
additional flux 6^ vanishes. The corresponding vector 
potential is of the form 



1 



(A8) 



where the overall prefactor has been chosen so that 
5B{R) = rj^o/nR^, i.e. the field strength at the surface 
of the cylinder is the same as it would be in the case of 
a uniform magnetic field. 

Inclusion of the above vector potential leads to the 
following additional term in the Hamiltonian, 



6h = CT3Si2ry 



i?2 



1 



(A9) 



Any cylindrically shaped surface with slight rotation- 
ally invariant deformations from ideal cylinder can be 
described by the function R = R{z) in which z is the 
distance from an arbitrarily chosen origin on the z axis 
(the axis of the deformed cylinder) and R{z) — i?o is the 
deviation from the ideal cylinder with radius Rq. One 
can use Eq.Q to find the Hamiltonian for the surface 
electrons of a TI wire with such a configuration if R{z) 
and its first derivative is continuous. In order to do so 
it is convenient to use a parameter, (3, which is defined 
as the angle between the normal vector n and the plane 
perpendicular to the axis of the wire. Note that we have 
assumed that the deformation from the ideal cylinder is 
such that the axis of the initial ideal cylinder always re- 
mains inside the wire and there is a rotational symmetry 
around that axis. This way we can write the normal unit 
vector in the cylindrical coordinate system 



n((/3, z) = cos (3{z)p{ip) -f sin /3(z)z 



(Bl) 



for the ideal cylinder /3(z) ~ (region A as shown in Fig. 
([5])). /3 gradually goes to ±7r/2 at the smooth ends of the 
wire. Using Eq. g we can obtain the continuum real 
space representation of the 2x2 Hamiltonian matrix for 
the surface states in the presence of the magnetic flux 



1 

hv 



ho = ^icff{z)+s^I)i{z)+Sp'D2{(p, z) + S3l)3{ip, z) (B2) 



in which s^p = —si sinip + S2 cos(/j and Sp = si cos(y9 
S2 sin if and the functions above are defined as 



2 \R{z) 



d(i{z] 
dz 



cos j3{z) 



(B3) 



Di(z)=zcos/3(z)5. + ^sin/?(z) (^-J-^ - ^ ) (B4) 



^2{^,z)=''-^^^{id^ + r,{z)) (B5) 



R{z) 
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cos/3(z) 



9sfJ-B 

hv 



Bo (B6) 



r]{z) is the fraction of the magnetic flux in the units of 
h/e that is enclosed by the radius R{z). 

The same unitary transformation that has been used 
in defining the spinors given in Eq. Q can transform 
isp,Sip) to (si,S2) and therefore it can make the new 
Hamiltonian matrix invariant under rotation (p ^ ip + 2TT 



'I 



,0 e'f , 



(B7) 



zero. Therefore in the region A and for 77 = 0.5 and I = 
we get 

hA{z)=iS2d, (B14) 

The region B Hamiltonian can be obtained similarly. 
Here we have (3 = 9 where 6 is the angle that 
parametrizes the quarter circle cross the section of the 
surface in region B. It varies from zero at the A/B 
boundary to 7r/2 at the B/C boundary. We thus have 



R{z) = Rq — ao{l — cos 9) 
z = gq sin 9 



(B15) 
(B16) 



One has to be careful with the partial derivative with 
respect to tp since it does not commute with the trans- 
formation matrix U {(p) . It is easy to check that it trans- 
forms as follows 



S3 
2 



(B8) 



This way the transformed Hamiltonian is invariant under 
rotations about the axis of the wire with the associated 
quantum number I. Therefore, similar to our treatment 
of the infinite wire, for each Z = 0, ±1, . . . we obtain 



1 r 

hv 



hi{z) = /iefF(z) + S2D(z) + simi/(z) + S3m2((z) (B9) 



The functions used in the above expression are defined 
as 



1 d/3 

r'^'^Tz 



(BIO) 



D{z) = -{d, cos (3 + cos pd,) (Bll) 



niu(z) = -^a+^-r;(z)) (B12) 



COS /3 1 

m2/(z) = -^^(/ + - - 7;(z)) + TOZccman (B13) 

One can recover the infinite cylinder Hamiltonian given 
in Eq. ([s]) by replacing R{z) with Rq and putting (3 to 



Assuming that ag ^ i?o and using 9 as our new coordi- 
nate we can obtain the Hamiltonian in this region for the 
case where half quantum magnetic flux is penetrating the 
interior of the wire 



hB{9) 



1 

2ao 



1 



is2de 



(B17) 



To obtain the Hamiltonian in the region C one has to 
start from Eq. ([T]) again since R{z) is not well deflned 
in this region. Using polar coordinates (p, (p) the Hamil- 
tonian takes a simple form and one can use the same 
unitary transformation given in Eq. ( |B8[ ) to put it in a 
rotationally invariant form. For the case where the half- 
quantum magnetic flux is uniformly distributed over the 
disk surface the I = Hamiltonian becomes 



hc{p) = -is2dp 



2p 



2) + S3TOZ0 



(B18) 



Now since we assume that the SC order parameter van- 
ishes in the C region the zero-energy solutions to the BdG 
Hamiltonian in this region are the same as solutions of 
the above Hamiltonian. Assuming furthermore that the 
Zeeman field is negligible, we obtain a first order differ- 
ential equation for each component of the spinor which 
is easy to solve and it leads to the solution of the form 
given in Eq. ( 40 1 



For regions A and B the presence of the SC gap leads 
to four coupled linear differential equations which can 
be decoupled by a linear unitary transformation. What 
remains is to match the solutions at the two boundaries 
A/B and B/C and the result of this straightforward but 
somewhat tedious procedure is presented in Eqs. (38p0|. 
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